31. Navier Stokes Equations

Find velocity \(u : \Omega \times [0,T] \rightarrow R^d\) and pressure \(p : \Omega \times [0,T] \rightarrow R\) such that

\[\begin{split} \begin{array}{ccccl} \frac{\partial u}{\partial t} - \nu \Delta u + u \nabla u & + & \nabla p & = & f \\ \operatorname{div} u & & & = & 0 \end{array} \end{split}\]

from ngsolve import *
from ngsolve.webgui import Draw

Schäfer-Turek benchmark geometry:

from netgen.occ import *

shape = Rectangle(2,0.41).Circle(0.2,0.2,0.05).Reverse().Face()
shape.edges.name="wall"
shape.edges.Min(X).name="inlet"
shape.edges.Max(X).name="outlet"
Draw (shape);

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.07)).Curve(3)
Draw (mesh);

Higher order Taylor-Hood element pairing:

V = VectorH1(mesh,order=3, dirichlet="wall|cyl|inlet")
Q = H1(mesh,order=2)
X = V*Q

u,p = X.TrialFunction()
v,q = X.TestFunction()

nu = 0.001  # viscosity
stokes = (nu*InnerProduct(grad(u), grad(v))+ \
    div(u)*q+div(v)*p - 1e-10*p*q)*dx

a = BilinearForm(stokes).Assemble()

# nothing here ...
f = LinearForm(X).Assemble()

# gridfunction for the solution
gfu = GridFunction(X)

parabolic inflow at inlet:

uin = CoefficientFunction( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))
Draw (gfu.components[0], mesh, "vel");

solve Stokes problem for initial conditions:

inv_stokes = a.mat.Inverse(X.FreeDofs())

res = f.vec - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res

Draw (gfu.components[0], mesh);

implicit/explicit time-stepping:

\[ \frac{u_{n+1}-u_n}{\tau} - \nu \Delta u_{n+1} + \nabla p_{n+1} = f - u_n \nabla u_n \]

and $\( \operatorname{div} u_{n+1} = 0 \)$

tau = 0.001 # timestep

mstar = BilinearForm(u*v*dx+tau*stokes).Assemble()
inv = mstar.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")

the non-linear convective term \(\int u \nabla u v\)

conv = BilinearForm(X, nonassemble=True)
conv += (Grad(u) * u) * v * dx

implicit Euler/explicit Euler splitting method:

t = 0; i = 0
tend = 5
gfut = GridFunction(V, multidim=0)
vel = gfu.components[0]
scene = Draw (gfu.components[0], mesh, min=0, max=2, autoscale=False)

with TaskManager():
    while t < tend:
        res = conv.Apply(gfu.vec) + a.mat*gfu.vec
        gfu.vec.data -= tau * inv * res    

        t = t + tau; i = i + 1
        if i%10 == 0: scene.Redraw()
        if i%50 == 0: 
            gfut.AddMultiDimComponent(vel.vec)
            print(f"t = {t}", end='\r')

t = 0.05000000000000004
t = 0.10000000000000007
t = 0.1500000000000001
t = 0.20000000000000015
t = 0.25000000000000017
t = 0.3000000000000002

t = 0.35000000000000026
t = 0.4000000000000003
t = 0.45000000000000034
t = 0.5000000000000003
t = 0.5500000000000004
t = 0.6000000000000004

t = 0.6500000000000005
t = 0.7000000000000005
t = 0.7500000000000006
t = 0.8000000000000006
t = 0.8500000000000006
t = 0.9000000000000007

t = 0.9500000000000007
t = 1.0000000000000007
t = 1.0499999999999952
t = 1.0999999999999897
t = 1.1499999999999841
t = 1.1999999999999786

t = 1.2499999999999731
t = 1.2999999999999676
t = 1.3499999999999621
t = 1.3999999999999566
t = 1.449999999999951
t = 1.4999999999999456

t = 1.54999999999994
t = 1.5999999999999346
t = 1.649999999999929
t = 1.6999999999999236
t = 1.749999999999918

t = 1.7999999999999126
t = 1.849999999999907
t = 1.8999999999999015
t = 1.949999999999896
t = 1.9999999999998905
t = 2.0499999999998852

t = 2.0999999999998797
t = 2.1499999999998742
t = 2.1999999999998687
t = 2.249999999999863
t = 2.2999999999998577

t = 2.349999999999852
t = 2.3999999999998467
t = 2.449999999999841
t = 2.4999999999998357
t = 2.54999999999983
t = 2.5999999999998247

t = 2.649999999999819
t = 2.6999999999998137
t = 2.749999999999808
t = 2.7999999999998026
t = 2.849999999999797
t = 2.8999999999997916

t = 2.949999999999786
t = 2.9999999999997806
t = 3.049999999999775
t = 3.0999999999997696
t = 3.149999999999764
t = 3.1999999999997586

t = 3.249999999999753
t = 3.2999999999997476
t = 3.349999999999742
t = 3.3999999999997366
t = 3.449999999999731
t = 3.4999999999997256

t = 3.54999999999972
t = 3.5999999999997145
t = 3.649999999999709
t = 3.6999999999997035
t = 3.749999999999698
t = 3.7999999999996925

t = 3.849999999999687
t = 3.8999999999996815
t = 3.949999999999676
t = 3.9999999999996705
t = 4.049999999999687
t = 4.099999999999704

t = 4.149999999999721
t = 4.199999999999737
t = 4.249999999999754
t = 4.299999999999771
t = 4.349999999999787
t = 4.399999999999804

t = 4.449999999999821
t = 4.4999999999998375
t = 4.549999999999854
t = 4.599999999999871
t = 4.649999999999888
t = 4.699999999999904

t = 4.749999999999921
t = 4.799999999999938
t = 4.849999999999954
t = 4.899999999999971
t = 4.949999999999988
t = 5.000000000000004

In the multidim - GridFunction gfut we have collected several time-steps, which can be animated now by the visualization:

Draw (gfut, mesh, interpolate_multidim=True, animate=True);